function gabormovie(write_avi_movie)

%for phase = -pi:pi/16:pi
    
%    imagesc(gaborfilter1((-2*pi+phase):pi/8:(2*pi+phase),10*pi,10*pi,pi/4)),colormap(gray)
%    drawnow
%end

fig=figure(100);
set(fig,'DoubleBuffer','on');

if write_avi_movie 
	mov = avifile('stimuli.avi')
end

phase = 0;
angle = 3*pi/4;
g=gaborfilter((-4*pi+phase):pi/3:(4*pi+phase),10*pi,10*pi,angle);
dc = .1;
for contrast = [-1:dc:1 1:(-dc):-1]
    
    imagesc(contrast*g,[-2 2])
    colormap(gray);
    drawnow;
    
    if write_avi_movie 
    	F = getframe(gcf);
        mov = addframe(mov,F);
        disp('out');
    end

    
end

if write_avi_movie
	mov = close(mov);
end

% Modified by Jeff McKinstry, 10-28-09
% Allows passing in an array of x values for the cross-section of the Gabor
% function.  This allows you to set the phase and the number of pixels in
% your output array.  
% 
% Sample usage:
% imagesc(gaborfilter1(-2*pi:pi/4:2*pi,pi,10*pi,pi/4)),colormap(gray)
%
%The Gabor filter is basically a Gaussian (with variances sx and sy along x and y-axes respectively)
%modulated by a sinusoid (with centre frequencies U and V along x and y-axes respectively) 
%described by the following equation
%%
%                            -1     x' ^     y'  ^             
%%% G(x,y,theta,f) =  exp ([----{(----) 2+(----) 2}])*cos(2*pi*f*x');
%                             2    sx'       sy'
%%% x' = x*cos(theta)+y*sin(theta);
%%% y' = y*cos(theta)-x*sin(theta);

%% Describtion :

%% one_d_fun : Input array of x values for the cross-section of the Gabor function; if n==length(one_d_fun) then G is nXn.
%% Sx & Sy : Variances along x and y-axes respectively
%% theta : The orientation of Gabor filter

%% G : The output filter as described above



%%  Author : Ahmad poursaberi  e-mail : a.poursaberi@ece.ut.ac.ir
%%          Faulty of Engineering, Electrical&Computer Department,Tehran
%%          University,Iran,June 2004

function G = gaborfilter(one_d_fun,Sx,Sy,theta);

% rotate
[x y] = meshgrid(one_d_fun);
xPrime = x * cos(theta) + y * sin(theta);
imagesc(xPrime)
yPrime = y * cos(theta) - x * sin(theta);
G = exp(-.5*((xPrime/Sx).^2+(yPrime/Sy).^2)).*cos(xPrime);

%G=double(real(G));

